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Abstract 

We study surface diffusion in the framework of a generalized Frenkel- 
Kontorova model with a nonconvex transverse degree of freedom. The model 
describes a lattice of atoms with a given concentration interacting by Morse- 
type forces, the lattice being subjected to a two-dimensional substrate poten- 
tial which is periodic in one direction and nonconvex (Morse) in the transverse 
direction. The results are used to describe the complicated exchange-mediated 
diffusion mechanism recently observed in MD simulations [J.E. Black and 
Zeng-Ju Tian, Phys. Rev. Lett. 71, 2445-2448 (1993)]. 
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I. INTRODUCTION 



Diffusion of atoms adsorbed on crystal surfaces is important in many processes such 
as surface reactions and crystal growth of artificially layered materialsBi. Besides, surface 
diffusion is of considerable intrinsic interest, because experimental results show a very rich 
and complicated behavior of diffusion coefficients, especially as a function of the atomic 
concentration. 

When the first adsorbed layer is complete, new incoming atoms start to fill the second 
adlayer. Usually the diffusion in the second layer follows the same laws as that in the first 
layer because the adatoms of the first monolayer play the role that substrate atoms played for 
the diffusion of the adatoms of the first layer, i.e. the first-layer adatoms create an external 
potential for the second-layer atoms. However, for some adsystems the situation may be 
more complicated owing to exchange of atoms between the first and second adlayers. Such 
an exchange was observed experimentally by Medvedev et all for the Li-W(112) and Li- 
Mo(112) adsystems, where the growth of the second layer results in the reconstruction of the 
underlying first adlayer. Moreover, recently Black and Tiani have observed a "complicated 
exchange- mediated diffusion mechanism" in a molecular dynamics experiment, where an 
isolated Cu adatom, which is diffusing on the Cu(100) surface, may enter the first substrate 
layer and create there a strain along a close-packed row. This localized excitation moves 
along the row for a distance of several lattice constants, and then the strain is relieved by an 
atom in the strained row popping out and returning to the surface. The simulation showed 
that this diffusion mechanism becomes important at high enough temperatures (T ~ 900 K 
for the Cu-Cu(lOO) adsystem). 

The effect of reconstructive crystal growth was studied theoretically in ref |6| within the 
framework of a generalized Frenkel-Kontorova (gFK) model with a nonconvex transverse 
degree of freedom. The model describes a chain of atoms interacting with a generalized 
Morse potential. The atoms are assumed to be mobile in two directions, one is along 
the surface (this is the direction along which the atoms can diffuse), and the other one 
is orthogonal to the surface. They are subjected to a two-dimensional substrate potential 
which is periodic along the chain (i.e., along the surface) and has the shape of a Morse 
potential in the transverse direction (orthogonal to the surface). The concentration of atoms 
is characterized by the dimensionless parameter 9 = N/M, the so-called coverage in surface 
physics, where N is the number of atoms and M is the number of minima of the external 
potential. This model may be considered as a "minimal" model which takes into account all 
main features of real adsorbed systems such as layers adsorbed on furrowed crystal surfaces. 
On the other hand, it is simple enough to be tractable analytically, at least in some aspects. 
The aim of the present paper is to show that the same models can provide a useful framework 
to study the role of the atomic exchange between the first and second adlayers in surface 
diffusion at coverages 6 > 1. 

In the study of surface reconstruction |6], it has already been shown that the formation 
of a metastable defect (kink) in which an atom of the second layer penetrates into the first 
layer is possible in some parameter range of the gFK model. As kinks of the first layer 
generally have a diffusion coefficient Dk which is much larger than the diffusion coefficient 
D act of thermally activated atoms of the second layer, it was speculated that the formation 
of such defects could have a large influence on the overall diffusion of atoms on crystal 
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surfaces. The present work confirms this conjecture and provides quantitative evaluations 
of the role of this exchange-solitonic mechanism on surface diffusion. We study a system 
of adsorbed atoms which has a complete first adlayer and an extra adatom that can be 
in two stable configurations. In the first one, the extra atom is in the second layer. This 
state corresponds to the minimum of the total potential energy (the gs-configuration). The 
second configuration with all adatoms in the first layer corresponds to a metastable state 
(the ms-configuration). In this ms-configuration the extra atom creates a localized solitonic 
excitation, the so-called kink, which usually has a very high mobility. Thus, even if the 
lifetime of the ms-configuration is small at a nonzero temperature, the ms-configuration 
may play the main role in surface diffusion, generating an unusual temperature dependence 
of the mass diffusion coefficient. 

The paper is organized as follows. The model is described in Section O. Calculations of 
quasi-adiabatic trajectories are described in Section |HI| , and the results of molecular dynam- 
ics simulation are presented in Section [TV[ Section [V] is devoted to theoretical estimations 
of different diffusion mechanisms. The last Section [VI] concludes the paper. 



II. MODEL 

We use the generalized Frenkel-Kontorova model introduced in ref []. The displacement 
of an atom is characterized by two variables: x describes its motion parallel to the surface 
and y describes its deviation orthogonal to the substrate. The potential perpendicular to 
the surface has a Morse shape 

V y (y)=e d (e'^-l) 2 , (1) 

which tends to the finite limit (known as the adsorption energy) when y — > oo. The 
parameter 7 determines the anharmonicity and it is related to the frequency u y of a single- 
atom vibration in the normal direction by the relation u 2 = 2 / y 2 e ( [/i r n, m being the atomic 
mass. It should be noticed that the function ([!]) is nonconvex, i.e. it has an inflexion point at 
V — Uinf = 7 _1 l n 2, which has important consequences on the properties of the system!. To 
model the substrate potential along the surface, we use the deformable periodic potential, 

= 1 (l + g ) 2 [l-cos(2Wa)] 
y ' 2 1 + s 2 - 2s cos(27rx/a s ) v ; 



as discussed in a previous workH. Here e s is the activation energy for diffusion of a single 
atom, a s is the period of the substrate potential along the chain, and the parameter s 
(|s| < 1) determines the shape of the potential. The frequency uo x of a single-atom vibration 
along the chain is connected to the shape parameter s by u 2 = Uq [(1 + s)/(l — s)] 2 with 
ojq = 27i 2 e s /ma 2 . In the following we use a system of units where a s = 2n, e s = 2 and 
m = 1, so that ujq — 1. 

The total potential energy of a single atom interacting with the substrate is written as 

V sub (x : y) = V x (x)e-i' y + V y (y). (3) 

The exponential factor in the first term of the right-hand side of Eq. (|3p takes into account 
the decrease of the influence of the surface corrugation as the atoms move away from the 
surface, so that V su b(x,y) — > Ed when y — » 00. 
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To model the interaction between adatoms, we use the generalized Morse potential 

Vmt[) a \(3-f3' 6 (3-(3' € /' [) 

where e a is the interatomic bonding energy of a molecule adsorbed on the surface, r a is the 
molecule's equilibrium distance, and the exponents (3 and (3' are related to the frequency oj a 
of interatomic vibration by the relation = e a /3/3'. 

We use the same model parameters as in the study of surface reconstruction [| Namely, 
considering the case of metal atoms adsorbed on a metal substrate such as, e.g., the Li- 
W(112) or the Li-Mo(112) adsystems, we take a s = 2.73 A which is the distance between 
the wells along a furrow on the W(112) surface, and e s ~ 0.1 eV, e& ~ 3 eV, cu x ^ uj y ~ 10 3 
cm -1 which are typical values for these systems. Returning to our system of units, we 
get e d = 60, 7 = 0.183, y in f = 3.80, u x = 1.5, u y = 2, 7' = 27 = 0.366, and s = 0.2. 
The interaction energy between two adsorbed metal atoms usually lies within an interval 
e a ~ 0.1 — 0.5 eVM, or e a ~ 2 — 10 in our system of units; we have chosen the value e a = 6. 
For the exponents (3 and (3' we take (3 = 1.9 and /3' = 0.19 (see a more detailed discussion of 
such a parameter choice in |6]). In the papers we had chosen for the interatomic equilibrium 
distance the value r a « 3.04 A (the interatomic distance in lithium metal), or r a = 7 in 
our system of units, because we had in mind the application of the model to the lithium 
film. However, we will use in the present work lower values of r a , r a = 6.3 and 6.4 in 
order to investigate the case when an extra atom can be inserted into the first adlayer and 
exist there in a metastable state. Although we have selected some of the parameters by 
comparison with a real system, we do not claim to describe quantitatively a concrete system 
of adatoms with a model which is still oversimplified. We are interested in the phenomenon 
of exchange-mediated diffusion, and this is why the parameter r a has been adjusted to allow 
such a phenomenon. 

In the present work we study the case of a fixed concentration of atoms. Therefore we 
impose periodic boundary conditions with a fixed number M of minima of the substrate 
potential as well as a fixed number N of adatoms (we have used M — 16 and iV = 17). 
The ground state configuration as well as the nearest metastable states are searched for 
with a standard molecular dynamics (MD) algorithm. Namely, we are starting from an 
appropriate initial configuration and allow the atoms to relax to a nearest minimum of the 
total potential energy of the system. Thus, the computer algorithm reduces to the solution 
of the equations of motion which follow from the potentials @ and (|J) with an artificially 
introduced viscous friction!!!. In computer simulations we can only include the interaction 
with a finite number of neighbors. This is achieved by introducing a cutoff distance r* 
(we have chosen r* = 5.5 a s ) and accounting only for the interactions between the atoms 
separated by distances lower than r* as usual in MD simulation. 

III. QUASI- ADIABATIC TRAJECTORIES 

Adiabatic characteristics of the system were studied in order to determine the potential 
energy surface of the model. First, we have determined the ground state with the MD al- 
gorithm with friction, in the same way as in the study of surface reconstruction []. Then, in 
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order to investigate mass diffusion, we have investigated the dependence of the potential en- 
ergy upon the position of the center of mass of the system. It is however difficult to impose a 
constraint on the center of mass. Therefore, instead of determining this adiabatic trajectory, 
we explore the multidimensional potential energy surface along "quasi-adiabatic" trajecto- 
ries defined by imposing appropriate constraints to a single atom. The most important 
characteristics of the motion, which are the positions and energies of the stationary points 
(the ground state, the metastable state, and the saddle points) are calculated correctly by 
this method, and we may expect to obtain the correct shape of the adiabatic trajectory at 
least qualitatively with this approach. 

Since we consider the possibility that the extra atom may exist in two stable states that 
differ by its distance to the substrate (atom in the second or in the first layer), a first analysis 
has been performed by constraining only the y coordinate yj of the extra atom. We displace 
it up and down in the y direction (perpendicular to the surface) by small steps. At each 
step we allow for the x coordinate of this atom and for both coordinates of all other atoms 
to adjust themselves to the new value of yj. For each relaxed state, the position Y of the 
center of mass of the system is calculated. Fig. [I] shows the dependence of the system energy 
upon Y. The ground state configuration corresponds to a Imposition of the kink around 2.8, 
while the metastable state corresponds to a much lower value of Y. These pictures allow us 
to calculate e ms , the difference in energies of the metastable and groundstate. The results 
are listed in Table Q. One can notice that e ms depends very much on r a ; the metastable 
state is much more likely to be relevant in the case (a) r a = 6.3 which has a lower e ms than 
in case (b) r a = 6.4. The two stable states are separated by a barrier Sbarrier- 

Around each value of Y corresponding to a stable state, we have then determined the 
potential energy as a function of the displacement along x of the coordinate Xj of the 
extra atom. In this CclSG Xj IS constrained to a given value but yj is allowed to relax to its 
equilibrium value as well as the x and y coordinates of all the other atoms. The X coordinate 
of the center of mass is again calculated for each relaxed state. The results are presented 
in Fig. for the case of a quasi-adiabatic trajectory starting from the ground state, and 
in Fig. for a trajectory starting from the metastable state. In this case, the structure of 
the metastable defect is such that the extra atom does not play a specific role and cannot 
be identified unambiguously. This has no consequence on the mass diffusion, but in the 
algorithm one has to chose the atom that will be constrained to a given Xj. The choice has 
been done in the same way as in refi. Figure |2] shows that, for the translation of the center 
of mass along x in the case of an extra atom in the second layer, there is no metastable 
state and only a barrier to overcome. The diffusion in the A-direction will therefore be an 
activated process and we call e act the difference between the unstable maximum state and the 
stable one. As attested by the results presented in Table [I], e act decreases when r a increases. 
More precise consequences for the physics will be explained in the following section. Figure 
|3] and Table | show that the barrier e pn for the A translation of the metastable state is 
extremely low compared to e act . 

In order to complete the picture of the potential energy surface, in a second series of 
calculations we artificially moved the same atom j, but now both coordinates Xj and yj 
were constrained, while the other atoms were allowed to relax to the minimum of the system 
energy. The atom j was moved to scan the x and y directions as in a TV sweep. The results 
allow to draw the full picture of the energy E of the system as a function of X and Y. They 
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are presented in Fig. f|a as a contour map and in Fig. |]b as a three-dimensional surface (we 
show the results only for r a = 6.3 , because the case r a = 6.4 looks qualitatively similar). 

^From the energy surface of Fig. |] and the quantitative results of Table [I], one can predict 
the general behavior of the system. Starting from the ground state in the second layer, the 
extra adatom may overcome the barrier e ac t and directly jump to the nearest neighboring 
site in the second layer, or it may overcome the barrier Sbarrier and form the metastable 
state in the first adlayer that will move practically without barriers. Therefore during the 
lifetime of the metastable state the local compression of the first layer may move for a long 
distance, and when the system returns back to the ground state an atom arises in the second 
layer in a site which may be far away from the initial position of the extra atom. Such an 
exchange-mediated diffusion should clearly be the favorable diffusion pathway in the case 
r a = 6.3 because, in this case, the activation barrier for the transition to the metastable 
state is lower than the activation energy in the second layer, starrier < £act- On the other 
hand, in the case r a = 6.4 we have the opposite inequality, Ebarrier > £ ac t, and direct jumps 
to a nearest site of the second layer should be the most probable diffusion path. However, 
even in this case, transitions to the metastable state may take place and, owing to the high 
kink mobility, these transitions may lead to a remarkable contribution to the total diffusion 
coefficient. This aspect is discussed in the following section. 

IV. MOLECULAR DYNAMICS SIMULATIONS. 

The predictions presented above are simply based on static configuration energies and 
collective dynamical effects could enter in a non trivial way to affect the diffusion. Thus 
these predictions must be checked with full molecular dynamics simulations of the system. 
Temperature effects are introduced through a standard Langevin approach. The hamiltonian 
equation of motion of each atom of mass m is completed by a friction term with a damping 
coefficient mi] and a random force with a 5 correlation function of factor ImriksT . Starting 
from the ground state configuration a first time interval t t h is allowed to reach the thermal 
equilibrium state. Then, during a time interval t run we save the dependencies X(t), Y(t), 
which define the position of the localized excitation and y max {t) which defines the maximum 
y-coordinate of the atoms at time t. This last quantity allows us to determine if the system 
is in the ground state or in the metastable state. 

A first step of the analysis has been devoted to the properties of the metastable at 
non zero temperature. In principle the position Y of the center of mass of the system 
should tell us whether the system of adatoms is in the ground state or the metastable state. 
However, in practice, measures of Y do not answer this question because the histogram for 
the probability P(Y) that Y takes a given value shows only one maximum corresponding to 
the ground state configuration. This does not mean that the metastable is never excited, but, 
because of the importance of the fluctuations it is not possible to distinguish the two stable 
states, characterized by the shift between the two layers of only one atom (among N=17 
atoms): the increase of Y is too small. The information on the state of the system can be 
deduced from y ma x because this quantity takes a large value if an atom is in the second layer 
(i.e. when the system is in the ground state). The histogram of the probability P(y m ax) 
does exhibit two well-defined maxima corresponding to the ground state and metastable 
state configurations up to temperatures above T = 2 (see Fig. |5|a). From this result, it is 
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possible to define an effective free energy of the system as 

F eff = -k B TlnP(y max ). (5) 
which is plotted in Fig. at different temperatures. The difference in free energy between 



the metastable state and the ground state (called e ms in Sec. 1I|), as well as the height 
^barrier of the barrier between the ground state and the metastable state are approximately 
equal to the values deduced from quasi-adiabatic calculations at T = 0. The relative vertical 
positions of the curves are not significant, because they depend on the reference in energy 
that we chose for each temperature. Note, however, that the free energy of the metastable 
state with respect to the ground state free energy, e ms (T), increases with temperature from 
the value 0.6 at T = 0.3 to 1.2 at T = 0.9. This effect may be understood by noticing 
that the variation of this difference in free energy is related to the entropy of both states as 
follows: 

derrJX) _ d[F eff (y^ ax )-F eff (y^ ax )] _ 

Qrp — Qrp — ^gs dms- {0) 

Since this variation is positive according to the numerical result, it means that the entropy 
of the ground state is higher than that of the metastable state. This results is reasonable 
because the atom isolated in the second layer has more space available to move than when 
it is in the first layer, and moreover its interaction energy with its neighbors and with the 
substrate (because of the exp(— ~fy) factor) is weaker. 

Integrating the peaks of the histogram of P(y max ), we can obtain the occupation numbers 
of the two states. The density of state for the ground state, p gs , is calculated by considering 
the states on the right of the minimum value between the two peaks. The density of state 
for the metastable state is then p ms = 1 — p gs . In Fig. ||b they are shown as functions of the 
temperature, and compared with the functions 

P 9 °J = (l + e-^^Y 1 and pffi = l-p£>, (7) 

which would be deduced from the quasi-adiabatic results without taking into account the 
variation of the energy of the metastable state with temperature. A significant deviation 
between the numerical values and the estimation provided by the quasi-adiabatic results 
shows up only in the high temperature range. 

The second set of studies has been devoted to a direct determination of the mass diffusion 
coefficient from the thermalized molecular dynamics simulation. The result is deduced from 
a series of k simulations at each temperature (k = 20), each one extending over the time 
interval t run . From simulation i we extract the diffusion coefficient using 

= {{ Xi t0 + t )-X M?} 

(where to + 1 is lower than t run ). Then we average over the k simulations. The results are 
presented in Figs. ^ and j|, and discussed in the following section. 

Besides the total diffusion coefficient, we calculate also the kink diffusion coefficient (i.e. 
the diffusion coefficient of the system in the metastable state). For this, we started from the 
metastable state and monitor the maximum y max coordinate of the atoms: when an atom 
escapes from the first adlayer, we stop the run. A reliable numerical estimation of the kink 
diffusion coefficient is only possible at low temperature when the metastable has a sufficient 
lifetime. The results are presented in Fig. P. 
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V. DISCUSSION 



Let us now come to the main point, an analytical estimate of the role of kinks in the pro- 
cess of surface diffusion. As shown above, the system may evolve from one GS configuration 
to an another GS configuration with the atom occupying a neighboring site in the second 
adlayer by two pathways, either a direct atomic jump over the barrier e ac t or through a 
two-stage mechanism involving the metastable state in the first adlayer. For such a two-way 
process the total diffusion coefficient can be calculated by the method described by Kutner 
and SosnowskaLl. The result is trivial: 

D = P gs D act + PmsDkink, (9) 

where the occupation numbers p gs and p ms were defined above. 

The hopping diffusion coefficient D act may be estimated with the help of the seminal 
result of the Kramers theory^ 



where u* (respectively a>&) is the frequency of the linearized oscillations of the system near the 
ground (resp. unstable) state. As shown on Fig. |^, the potential energy for the motion over 
the barrier is well approximated by a sinusoidal function V(x) = e act [l + cos(27rx/a s )] /2, so 




that uj* = Ub = ye act /2m. Finally, as a s = 2ir and m = 1 with our units, we get 

2,exp(-||), (11) 

To estimate the diffusion coefficient Dk in k of the metastable state, recall that the atom 
inserted into the chain creates a local distortion that can be considered as a kink in the 
Frenkel-Kontorova model. In the continuum approximation the equation of motion reduces 
to the exactly integrable sine-Gordon (SG) equation. This limit, in which a defect would 
propagate freely as a soliton in the system, can be used as a starting point for a perturbative 
analysis of the influence of the discreteness of the lattice. When the continuous translational 
invariance is broken, one finds that the kink experiences a periodic potential, with the 
periodicity of the lattice, which is known as the Peierls-Nabarro barrier in dislocation theory. 
A precise evaluation of the PN barrier turns out to involve subtle effects, particularly when 
discreteness effects are weakEl, but a recent analysis has derived accurate results in this 
case". The basic results is that the shape of the kink can be obtained by solving a Klein- 
Gordon equation obtained by replacing the actual substrate potential V(x) by an effective 
potential 



V eff (x) = V(x) 



1 



2Ag 



V'{x) 



2 



(12) 



where g = a^l / /^(a s )/(27r 2 e s ) = V^ t (a s ) is the elastic constant of the atomic chain. This 
expression allows us to derive the value of the mass of the kink which is given by the general 
expression for a Klein Gordon model, 



S 



m 



m kink = ^—j ^2V eff (x)dx . (13) 
a sy 9 "'o 

This formula gives the usual expression rrikink = 8 m /(a%-Jg) m the sine-Gordon case if, 
instead of the effective potential, we use the actual potential V(x) = [1 — cos(27nr/a s )]. For 
the values of g which are obtained with our parameters (see Table H), the correction to the 
SG value due to discreteness is small, giving for instance rrikink = 7.95 mj (a 2 s y/g) for the case 
r a = 6.3. The proper treatment of discreteness is much more important for the evaluation 
of the amplitude of the Peierls potential which is obtained a 

e pn « 712.26 e s g exp(-7r 2 v ^) . (14) 

Table |TT| gives the values of g, rrikink and e pn which result from the interaction potential (f|) 
with r a = 6.3 and r a = 6.4. It shows that the amplitude of the PN potential is negligible 
with respect to the barrier e act corresponding to the translation of the extra atom in the 



second layer. This is in agreement of the numerical results of Sect. |TTT| . Moreover, this value 
Epn <C ksT at the temperature that we consider, shows that the kink motion is not thermally 
activated. Its diffusion coefficient can be derived from the formula for a free diffusion 

Dkink = , (15) 

TTlkink Vkink 

where rjkink is an effective viscous friction for kink motion. At T = r]ki n k is simply given 
by rjkink ~ Vj where r\ is the viscous friction for the motion of an isolated adatom due to 
energy exchange with the substrate. The value of r] is usually^ about r] ~ uo/10. But at 
nonzero temperature rj kin k depends on T because of the coupling of the lattice phonons with 
the highly nonlinear core of the kink. The simplest form for this dependence is@ 

Vkink ~ V + aT, (16) 

where a is a coefficient which cannot be obtained analytically but could be obtained exper- 
imentally as in the case of copper0. 

In our analysis, we treat rrikink and a as adjustable parameters. They are chosen by 
fitting with Eq. (|15"D the temperature evolution of D^nk determined by the MD simulations 
with r] = 0.1, r a = 6.3. The parameters rrikink = 0.0886 and a = 0.155 obtained in this 
particular case are also valid for the case rj = 0.3 as shown from Fig. ^. Moreover, the 
expression flT5] ) with the same parameters describes the simulation results for the diffusion 
in the first adlayer for the case r a = 6.4 as well. This rather good agreement between 
the numerical results and the analytical estimation of D^nk shows that, although the SG 
description may seem rather crude for the generalized FK model, it provides a good basis 
for analysis. This is confirmed by the comparison between the fitted value of rrikink and the 
theoretical value given by Eq. (|13|): the fitted value is smaller than the theoretical one, but 
the order of magnitude of rrikink is however correctly given by the discrete SG calculation. 

Having obtained analytical estimates for the two diffusion coefficients D act and Dkink, 
we are now in a position to estimate the role of the kink diffusion in the general process 
of atomic diffusion in the adsorbed layer. The role of both mechanism is well illustrated 
on Fig. 0. First one can observe that the temperature dependence of the total diffusion 
coefficient D deduced from the MD simulations is well reproduced by Eq. (M). It should 
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be noticed that, at this level the fit is obtained without adjustable parameters since our 
parameters have been deduced from the study of D kink alone and the density of states p gs 
and p ms are the measured quantities. The good agreement points out that the diffusion 
mechanism that we propose, involving two different processes, provides a correct description 
of diffusion as it can be observed in MD simulations of the generalized FK model. It is also 
important to notice that both mechanisms are essential to reproduce the numerical results. 
If one would consider only one of the two processes, i.e. assume that p gs = 1 or p ms = 1 , the 
diffusion coefficient would evolve versus T along the dash-dotted or dotted lines of Fig. [7]. 
The theoretical value of D would disagree completely with the MD results. 

Fig. |8] summarizes the results by providing a comparison between the temperature de- 
pendence of D observed by MD and the theoretical value of Eq. (|9]) for r a = 6.3 and r a = 6.4 
and two values of the damping coefficient r\. This figure shows that all the numerical results 
are well described by the two-process diffusion with only two adjustable parameters a and 
m kink- The results show also that, even in the case r a = 6.4, for which the direct jumps to 
a nearest adsite in the second layer has to overcome a much lower barrier than the barrier 
for the transition to the metastable state, the solitonic-exchange mechanism involving the 
metastable state still brings a crucial contribution to the total diffusion coefficient D. This 
is due to the extremely high mobility of the metastable state. 

VI. CONCLUSION 

Following our results, one can distinguish three different diffusion mechanisms for atoms 
adsorbed on a crystalline surface. 

The first one is the conventional diffusion (sometimes called hopping diffusion) which 
involves only one atomic layer. In can be an individual process when an adatom directly 
jumps from one adsorption site to a nearest adsite. The diffusion coefficient in this case 
follows the Arrhenius law D(T) = Aexp (—E ac t/kBT) with the preexponential factor being 
independent (or weakly dependent) on T. At high coverages 9 ;$ 1 (as well as for 6> < 2, 
etc) owing to the interaction between the adatoms, the activation energy e act decreases 
with respect to the height of the original substrate potential e s . For a strong interatomic 
interaction the motion takes a collective (concerted) character and may be described with the 
help of the kink concept^™. In this case the activation energy e act for the chemical diffusion 
becomes the Peierls-Nabarro barrier e pn which is significantly lower than the individual 
activation energy. Therefore the process is generally no longer thermally activated and 
D[T) is approximately given by the free diffusion formula D(T) « fc B T /rrikinkVkink- 

While for the conventional diffusion the underlying adatoms in the first adlayer play 
a passive role creating the effective external potential for the diffusion of atoms in the 
second adlayer, for the exchange diffusion mechanism this role becomes active. Two types 
of mechanisms can be considered. The "conventional" (or "one-step") exchange diffusion 
mechanism has been known for a long time. It occurs when an adatom A from the second 
layer "pushes out" an adatom B from its regular position in the first adlayer and occupies 
the free site so created. The new configuration has only a single adatom (B) in the second 
adlayer. The intermediate configuration with both adatoms A and B inserted into the first 
adlayer is unstable. It corresponds to a saddle point in the potential energy surface. The 
elementary diffusion step takes a short time ~ 10~ 13 sec. The diffusion coefficient for the 
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one-step exchange mechanism should follow the same Arrhenius law as for the conventional 
hopping diffusion with the activation energy corresponding to the saddle state. Exchange 
diffusion may be important for coverages 9 > 1, when the first adlayer is complete and the 
second one starts to grow. 

The new diffusion mechanism studied in the present work, is a two-step exchange diffusion 
mechanism. The main difference with one-step exchange diffusion is that the configuration 
with the adatom inserted into the first adlayer, now corresponds to a metastable state in- 
stead of the unstable (saddle) state. Because this metastable state corresponds to a kink 
configuration which is characterized by a very high mobility, the mechanism may be called 
as the "exchange-solitonic" mechanism of surface diffusion. For this mechanism the diffusion 
follows an Arrhenius law too, but now the activation energy corresponds to the difference 
in energies between the metastable state and the ground state (and not to the barrier for 
the transition from the second adlayer to the first one), while the preexponent factor is 
determined by the kink diffusion coefficient and essentially depends on temperature. The 
mass diffusion coefficient in this situation is described by the expression (|) where D act and 
D kink are given by Eqs. ([11]) and (|1|). 



We considered above the situation where an atom is diffusing in the second adlayer over 
a filled first adlayer. It is clear that the same situation can occur for adatoms of the first 
layer when the adatoms are the same as the atoms of the substrate. Indeed, in this case the 
top-layer substrate atoms play the same role as that of the first-layer adatoms in the former 
case. In particular the exchange diffusion mechanism was observed for the first time when 
De Lorenzi and Jacuccfl investigated by the MD method the self-diffusion of Na atoms 
adsorbed on the surface of a Na metal crystal. They found that the exchange mechanism is 
responsible for the diffusion across the rows on the furrowed crystal surface. This conclusion 
was later confirmed by the field ion microscope technique for the diffusion of atoms adsorbed 
on the transition metal surfaces^!. 

It seems possible and would be interesting to find and investigate the exchange-solitonic 
mechanism of surface diffusion experimentally with the field-ion microscope or STM tech- 
nique. It may be predicted that such a diffusion is to be expected for coverages 6 ~ a s /r a . 
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TABLES 



TABLE I. Parameters of quasi-adiabatic trajectories for two values of the equilibrium distance 
r a of the interatomic potential. e ms is the difference in energies of the metastable and ground state 
configurations, Ebarrier is the activation barrier for the transition from the gs-configuration to the 
metastable state, e ac t is the activation energy for the hopping diffusion, and e pn is the barrier for 
kink motion parallel to the surface. 



parameter 


r a = 6.3 


r a = 6.4 




0.755 


3.014 


^■barrier 


1.752 


3.476 


£act 


2.364 


2.055 




4.18 • 10- 


6 6.39 • 10- 6 


TABLE II. 


Parameters corresponding 


to the kink (metastable state): g is the elastic constant, 


fnkink is the kink mass, and £ pn (theory) is the amplitude of the Peierls-Nabarro barrier. 


parameter 


r a = 6.3 


r a = 6.4 


9 


2.243 


2.759 


™kink 


0.1344 


0.1214 


£pn(theory) 


1.22 • 10" 


3 2.98 • 10~ 4 
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FIGURES 



FIG. 1. Change of the energy AE with respect to the ground state as function of Y for a 
displacement along the quasi-adiabatic trajectory perpendicular to the surface for (a) r a = 6.3 and 
(b) r a = 6.4. 

FIG. 2. Change of the energy AE with respect to the ground state as function of X for a 
quasi-adiabatic displacement parallel to the surface starting from the ground state for (a) r a = 6.3 
and (b) r a = 6.4. 

FIG. 3. Change of the energy AE as function of X for a quasi-adiabatic displacement parallel 
to the surface starting from the metastable state for (a) r a = 6.3 and (b) r a = 6.4. 

FIG. 4. Potential energy surface for the case r a = 6.3. The X-Y contour map is plotted in 
figure (a) while figure (b) presents a three-dimensional plot. 

FIG. 5. Simulation results for r a = 6.3 and rj = 0.1. (a) Probability of the maximum of 
the y position for different temperatures. The solid line corresponds to T = 0.3, the dotted 
line corresponds to T = 0.5, the dashed line corresponds to T = 0.7, and the dot-dashed line 
corresponds to T = 0.9. (b) Density of states versus temperature. The solid line and stars describe 
the simulation results, whereas the dashed line corresponds to the estimation @. 

FIG. 6. The effective free energy (|5j) as function of the maximum of the y atomic positions as 
follows from the results of Fig. |^. The solid line corresponds to T = 0.3, the dotted line corresponds 
to T = 0.5, the dashed line corresponds to T = 0.7, and the dot-dashed line corresponds to T = 0.9. 

FIG. 7. Diffusion coefficient versus temperature for r a = 6.3 and ry = 0.1. The diamonds 
correspond to the simulation results. The error bars are computed numerically with 20 simulations. 
The dotted curve corresponds to Dkink and the dash-dotted curve corresponds to D act . The stars 
and the solid curve describe the estimation result D = p gs D act + p ms Dki n k-, where p gs and p ms were 
taken from the simulation results of Fig. ||b. 

FIG. 8. Comparison between the MD results and the theoretical values of the diffusion coeffi- 
cient in four cases. The symbols correspond to the simulation results, whereas the solid curves, to 
the theoretical ones. The diamonds describe the case r a = 6.3 and rj = 0.1, the triangles, the case 
r a = 6.3 and rj = 0.3, the squares correspond to the case r a = 6.4 and rj = 0.1, and the stars, to 
the case r a = 6.4 and rj = 0.3. 
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FIG. 9. The kink diffusion coefficient versus temperature for r a = 6.3. The symbols correspond 
to the simulation results and the lines describe the estimated results of Eq. (15). The stars and 
the dotted curve correspond to rj = 0.1, the diamonds and the dashed curve, to rj = 0.3. The 
parameters m^nk and a were obtained by fitting the expression ( Jl5| ) to the simulation data for the 
case rj = 0.1, and then the same values were used for ry = 0.3. 
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